library(data.table)
library(maps)
dir = '/Users/coryzigler/Dropbox/Research/HEI Accountability Grant/Data/'
realdatdir = '/Users/coryzigler/Data/Medicare/Linked AQS-Medicare/HEI PM10 Accountability/'

## unique site identifier = Monitor (FIPS-SiteID)
## unique monitor identifier = (Monitor, POC)

########################################
###      Medicare + AQS Data        ####
########################################
# These files were created with  Google Drive/ARP/Arepa/Link_PM10.R
#### If using fake data
## dat = as.data.frame(fread(paste(dir, 'Merge_6mile.csv', sep="")))  

#### If using real data
## dat = as.data.frame(fread(paste(realdatdir, 'Merge_6mile.csv', sep="")))

dat$FIPS = substr(dat$Monitor, 1, 5)


########################################
###  Attainment/Nonattainment Data  ####
########################################
attaindat=read.csv(paste(dir, 'EPA Nonattainment Table.csv', sep=''))
attaindat$FIPS=paste(formatC(attaindat$fips_state, width=2, flag="0"), formatC(attaindat$fips_cnty, width=3, flag="0"), sep="")
attaindat=subset(attaindat, pollutant=='PM-10')

### OPTIONS FOR DEFINING CROSS-SECTIONAL NONATTAINENT STATUS
##   0= attainment, 1=nonattainment if ever designated for PM10
#attaindat$a=as.numeric(!is.na(attaindat$pw_1990)) #attaindat$a = ifelse(attaindat$pw_1990 %in% c("P", "W"), 1, 0) #0= attainment, 1=nonattainment if designated in 1990

##   0 = attainment, 1 = nonattainment if diagnosed between certain years
whichyrs = 1991:1995
attainyrs = as.vector(attaindat[, paste("pw_", whichyrs, sep="")])
attaindat$a = ifelse(apply( attainyrs =="P" | attainyrs =="W", 1, sum)>0, 1, 0) #1=attainment in 1990-1995, 0=else

attaindat=attaindat[, (names(attaindat) %in% c("FIPS", "a"))]



########################################
###        Census Data               ###
########################################
censdat=read.table(paste(dir,'Cory-Census-2000-04-11-2011', sep=""), sep=",", header=TRUE)
names(censdat)[names(censdat)=="Master_county"]="FIPS"
names(censdat)[names(censdat)=="white_rate"]="White_cens"
names(censdat)[names(censdat)=="black_rate"]="Black_cens"
names(censdat)[names(censdat)=="Other_rate"]="Other_cens"
names(censdat)[names(censdat)=="Hispanic_rate"]="Hispanic_cens"
names(censdat)[names(censdat)=="Female_rate"]="Female_cens"
censdat$FIPS = as.character(formatC(censdat$FIPS, format="d", width=5, flag="0"))



########################################
###         Weather Data             ###
########################################
weather=read.csv(paste(dir,'Annual-weather-11-05-2012.csv', sep=""), sep=",", header=TRUE)
weather=subset(weather, year==1990)
weather$FIPS = as.character(formatC(weather$FIPS, format="d", width=5, flag="0"))




########################################
###            Merges                ###
########################################


d1 = merge(dat, attaindat, by="FIPS", all.x=TRUE)
d1$a[is.na(d1$a)] = 0

d2 = merge(d1, censdat, by="FIPS", all.x=TRUE)

d3 = merge(d2, weather, by = "FIPS", all.x=TRUE)



########################################
###            Output                ###
########################################
## outdir = realdatdir
## outdir = '/Users/coryzigler/Dropbox/Research/HEI Accountability Grant/PM10 Accountability/Data/'
alldat = d3
save(alldat, file = paste(outdir, 'alldata.RData', sep="" ))






